{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model selection using lmfit and emcee "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`lmfit.emcee` can be used to obtain the posterior probability distribution of parameters, given a set of experimental data. This notebook shows how it can be used for Bayesian model selection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import lmfit\n",
    "import matplotlib.pyplot as plt\n",
    "import corner\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# calculate a Gaussian\n",
    "def gauss(x, a_max, loc, sd):\n",
    "    return a_max * np.exp(-((x - loc) / sd)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Generate some data\n",
    "x = np.linspace(3, 7, 250)\n",
    "np.random.seed(0)\n",
    "y = 4 + 10 * x + gauss(x, 200, 5, 0.5) + gauss(x, 60, 5.8, 0.2)\n",
    "dy = np.sqrt(y)\n",
    "y += dy * np.random.randn(np.size(y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Container object of 3 artists>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl4VdX1/t8VQmROCJAQCPOMMgpI0WIqCk4FtUWxTnVs\na+2gtXVshTpUbbWTVavUVv22VutI/YmKSgRFQBAQCCggQwYykAQIhJn1++O9p3fInZLcOevzPPe5\n5+wz3HVPbt6zztprry2qCsMwDCN1SYu3AYZhGEZ0MaE3DMNIcUzoDcMwUhwTesMwjBTHhN4wDCPF\nMaE3DMNIcYIKvYi0EZFlIrJaRIpE5Deu9mwRWSAiX4rIuyKS5XHMHSKySUQ2isjUaH8BwzAMIzgS\nKo9eRNqpar2IpAP4CMCtAKYD2KWqD4vIbQA6q+rtIjIcwL8AjAfQE8B7AAar6vGofgvDMAwjICFD\nN6pa71rMANAKQC0o9M+62p8FcIFreQaAF1T1iKpuA7AZwIRIGmwYhmE0jpBCLyJpIrIaQAWAhaq6\nHkCuqla4dqkAkOta7gGgxOPwEtCzNwzDMOJEeqgdXGGX0SKSCeAdEfmGz3YVkWDxH6uxYBiGEUdC\nCr2Dqu4Rkf8H4GQAFSLSXVXLRSQPQKVrt1IAvTwOy3e1eRHixmAYhmEEQFWlsceEyrrp6mTUiEhb\nAGcBWAVgHoCrXLtdBeB11/I8ALNEJENE+gEYBGB5AGMT/nXPPffE3YZUsTMZbDQ7zc5EfzWVUB59\nHoBnRSQNvCk8r6rvi8gqAC+JyLUAtgG42CXeRSLyEoAiAEcB3KjNsc4wDMNoNkGFXlXXAhjrp70G\nwJkBjnkAwAMRsc4wDMNoNjYyNggFBQXxNiEsksHOZLARMDsjjdmZGIQcMBWVDxWxiI5hGEYjERFo\npDtjDcMwjOTHhN4wDCPFMaE3DMNIcUzoDcMwUhwTesMwjBTHhN4wDCPFMaE3jABMnw7U1sbbCsNo\nPib0hhGAhQuB7dvjbYVhNB8TesPww/HjwL59QEVF6H0NI9ExoTcMP+zfz/fy8vjaYRiRwITeMPxQ\nV8d38+iNVMCE3jD8YEJvpBIm9EaLZ+1a4PbbvdscobfQjZEKmNAbLZ7PPgMWLfJuq6sDRMyjN1ID\nE3qjxVNaCuzd691WVwfk55tHb6QGJvRGi6esDNizx7utrg4YNMg8eiM1MKE3Wjylpf6Fvm9fYPdu\n4MiRuJhlGBHDhN5o8ZSWUtiPH3e31dUBWVlAly5AVVX8bDOMSGBCb7R4ysr47mTaOMsdOwK5uRa+\nMZIfE3qjRXPsGFBZSUH3DN+Y0BuphAm90aKprASysxmi8cy82bePQt+9u2XeGMmPCb3RoqiqAg4d\ncq+XlgI9egCZmf49+r59gQ0bYm6mYUQUE3qjRXHFFcBTT7nXS0uBnj2BTp38C/355wNvvBF7Ow0j\nkpjQGy2GI0eAjz4C3nrL3VZW5vboPUM3jtCPG8cwjnn1RjJjQm+0GFatAnJyKPYHDrDN8egDhW7S\n0oALLgBeey0+NhtGJDChN1oMixYB550HjB4NfPgh2xyPvlMn/x49wCkFPZ8CDCPZCCr0ItJLRBaK\nyHoRWSciP3a1zxaREhFZ5Xqd43HMHSKySUQ2isjUaH8Bw/Dk8ceBrVv9b/vwQ2DyZODss4H589m2\naxfQrVtgjx4ABg4ESkq4vHUrUFwcPfsNIxqE8uiPALhZVU8EMBHAD0VkGAAF8KiqjnG95gOAiAwH\ncAmA4QDOBvC4iNhTgxEzXnyRIRpfVBmymTwZGD8eKCpi+65dQNeu/oW+QwcuO7n0qsAf/wg8+WT0\nv4dhRJKgIqyq5aq62rW8D8AGAD1dm8XPITMAvKCqR1R1G4DNACZEzlzDCE59PVBb616/4w7g00+Z\nL5+eTtHu1s1d1mDXLubQe4Zujh5lCmb79lxv356x+n37gJ07rSSCkXykh7ujiPQFMAbAUgCnAviR\niFwJYAWAn6nqbgA9XNsdSuC+MRhG1KmvZyEygGL92GOMwasCffqwvVs3Cj8AVFd7e/S//jUHULVv\nz3r0Do5XX1YGHD4c2+9kGM0lrLCKiHQA8DKAn7g8+ycA9AMwGsBOAI8EOVyba6RhhIun0C9cSC/8\nq6+A7du9hX7XLnrue/YAnTu7hf7NN4GHHnLH5x0cod+5032TMIxkIaRHLyKtAbwC4P9U9XUAUNVK\nj+1zAfzXtVoKoJfH4fmutgbMnj37f8sFBQUoKChonOWG4QdPoZ83D5gwgR2oPXu6hf6EE4C2bdme\nlQW0asXQza5dwJYtFP3MTO/z5uS4hT7Nep2MGFFYWIjCwsJmnyeo0IuIAPgbgCJV/YNHe56q7nSt\nXghgrWt5HoB/icijYMhmEIDl/s7tKfSGESkcoVel0P/2t8BvfgP06sXsGYecHA6C6tKF65mZwPr1\nwIknAt/8JrBggfd5c3OBzZuBgwctRm/EDl8neM6cOU06TyiP/lQAlwP4XEScXIY7AVwqIqPBsMxW\nAN8DAFUtEpGXABQBOArgRlW10I0RE1TdQl9dDezfD5xzDnDDDaxZM2WKe99u3Zh507Ur1zMzefy4\nccCNN3J/T3JzgTVrgH79GAY6fBjIyIjVNzOM5hEq6+YjVU1T1dGeqZSqeqWqjlTVUap6gapWeBzz\ngKoOVNWhqvpO9L+CYZDDhzl5yO7dDLN0787QTOvWwIoV7tAN4PboHaHv1Inv48dzrtgbbvA+d04O\n0zZ79uQxu3a5t+3fDzz6aHS/m2E0B4s2GilDfT3fHaHPzeV6//6MrXsKfbdu3qGbNm3ooY8b5//c\nubnAF18AeXneWTsAsG4dcOut7glMDCPRMKE3Uob6eqZE+gp9v34c/NS5s3tfX48eAH71K2DUKP/n\nzsnh00JeHperqpiff+wYR82qAq+8Er3vZhjNwYTeSGpqaymyAIU+J4dtvh59nz7eefHdujH10lPo\n77qLGTn+cM7l6dGffTawfDkLo/XvD/znP5H/foYRCUzojaTm3HPdGTKO0B88SPHNyWG7I/SedOvG\nd0+hD4Yj9D168NiPPwZqapijX1ICXHUVsHYtQ0SGkWiY0BtJy/HjFNdly7heX88RrVlZwJdfusV5\n5kymWHri3AScGH0osrJYQsEJ3cybx3ZH6AcMAMaOZbzeMBINE3ojadmxgxkvK1Zwvb4eaNeOovzF\nF26hz84GRo70PraxHn1aGrN4evbksaWlPKcj9Pn5zNW3ypZGImJCbyQtRUUcBLVyJdc9hX7LFrfQ\n+8Px6MMVegAoLASGDHEfe9llJvRGcmBCbyQtRUWM0R84AJSXewv90aPBhd4R+HBDNwDDMyL06Nu0\nAS68kKNlnVmq8vNN6I3ExITeSFqKiliy4OST6dV7Cj3g9rz9kZEB/PWvDOs0lpNOAu67j528VVUs\ngNamjXn0RuJiQm8kLUVFwPDhFPoVK7yFvkMHLgfjhhuaVqAsKwv42c9YDK1PH3rygAm9kbiEXY/e\nMBKFVas4LWBRETBsGMMnCxdS2J1XsLBNJBkwgNk4gFvoVb1z9g0j3pjQG0nH/Pmczq9/f8bYc3I4\ngMnx6DMyYif0/fszzRNwlzbes4di7zkS1zDiiYVujKSjvh64/npg9WquO1MDeoZuYiX0M2awrDFA\nL75XL+C114DRo2Pz+YYRDubRG0nH/v3enaieHn1+PssR+5YZjhbTpnmv9+oF3Hsvs4AChXBUGXqy\nuXaMWGEevZF07N/vnrgbcHv0+/fTox86FDjvvPjY1qsX69UfP860T39UVQFnnukO+RhGtDGhN5IO\nJ0Tj0K4ds2cqK0Nn2kSb3r1Z7Cwnh7Vw/FFezqqXdXWxtc1ouZjQG0mHr0cPUFi3bYu/0N94I/DU\nU+yIra31v49T+CzQjcAwIo0JvZF0+BP6bt0YMom30HftylGywYS+vJzv1dWxs8to2ZjQG0mHE4v3\nJCeH9eV9bwDxIjs7sMduHr0Ra0zojaTDKUfsiVONMt4evUM4Hr0JvRErTOiNpCNQjB5IHKEP5tGX\nl3NwlQm9EStM6I2kw1/oJpk8+p07WYzNYvRGrDChN5IOf6GbRPbo6+qA005zz21bXk6hN4/eiBUm\n9EbC8/bbwOuvu9cDZd0AiSP0nh59WRnnmN2+nes7d7Lqpgm9EStM6I2EZ8kS4P33uXzkCEeUZmR4\n79OtG8sNnHBC7O3zh6fQ79rF9xUreJM6fJjF0Cx0Y8QKE3oj4amvd4umE5/3rSGTk+O/PV54hm4c\noV+5kmGbvDxW3TSP3ogVJvRGwuMp9P7i8wAHKT35ZGztCoavR5+bS4++vJyTjAfLyjGMSGNCbyQ8\nBw64RdFffB7gbE+XXx5bu4Lh69FPnUqPvqyMHr0JvRFLTOiNhMdf6CbRycriBCTHj1PoTzqJdt90\nE3DOOW6P38nEMYxoElToRaSXiCwUkfUisk5EfuxqzxaRBSLypYi8KyJZHsfcISKbRGSjiEyN9hcw\nUp/6erf3Gyh0k2ikp1PY9+6l0HftCvz0pyx4dv317Exu25bbDSPahJp45AiAm1V1tYh0ALBSRBYA\nuBrAAlV9WERuA3A7gNtFZDiASwAMB9ATwHsiMlhVrfK20WQOHHB7v4FCN4lIdjbtdoT+mmsabq+p\ncU9BaBjRIqhHr6rlqrratbwPwAZQwKcDeNa127MALnAtzwDwgqoeUdVtADYDmBAFu40WRH09cPQo\nRT5ZQjeAOzzjCL0v2dmWYmnEhrBj9CLSF8AYAMsA5KpqhWtTBQBnhs4eAEo8DisBbwyG0WTq6/le\nU5NcHn1+PgdJmdAb8SasOWNdYZtXAPxEVevEI1lZVVVEgnUp+d02e/bs/y0XFBSgwCbQNAJw4AAH\nQtXWJk+MHmAH7Lp1FHN/Qt+5M7B7d+ztMpKHwsJCFBYWNvs8IYVeRFqDIv+8qjoD0StEpLuqlotI\nHoBKV3spgF4eh+e72hrgKfSGEYz6eqBHDwp9MoVuRowAXn6ZtW6yshpuz8xkZo5hBMLXCZ4zZ06T\nzhMq60YA/A1Akar+wWPTPABXuZavAvC6R/ssEckQkX4ABgFY3iTLDMNFfT0HRCVb6Oakk4DFi+m5\np/n5TzOhN2JFKI/+VACXA/hcRFa52u4A8CCAl0TkWgDbAFwMAKpaJCIvASgCcBTAjaqWKWw0jwMH\nKPTJFroZOpShmUGD/G83oTdiRVChV9WPENjrPzPAMQ8AeKCZdhkGAA44OniQoRvHo/cX705EMjKA\nwYNZ18YfWVnApk2xtclomdjIWCOhOXiQHbFOTnoyxegBxukD3ZgyM92dsRUV/vcxjEhgQm8kNAcO\nUNg9hT5ZQjcAMHKku1a+L07o5tAhoF8/4Nix2NpmtBzCSq80jHhRX0+h79yZoZtDh5JL6G+80T0O\nwBdH6KuqeEOrrwc6doytfUbLwITeSFiczte2bd0evWpyhW6ysvynVgLeQg8wDdOE3ogGFroxEpKN\nG4HJk92hm86dWeJ3587k8uiD4Qh9pWsUyr598bXHSF1M6I2EZNcuCrsTusnLY4bKsGHA6NHxti4y\nOKWMHY/ehN6IFha6MRKSujqGavbuZeimVy92xKan0C+2Uydvj76uLr72GKmLefRGQlJXx3h8aak7\nJp9KIg8ArVsDbdoA27Zx3Tx6I1qY0BsJQW0tMHeue93xbouLk6vztbFkZgKbN3PZPHojWpjQGwnB\nmjXAHzyqKTkzL23fztBNquIIfXa2efRG9DChNxKCXbu8p9VzvNsdO1Lbo8/KArZu5YApX49+zRp3\n/N4wmoMJvZEQVFd7C11LCt0cPQr079/Qo//1r4E33oiPXUZqYUJvJASO0Du1TuvqKII7dqR+6Aag\nR+8r9KWlFrc3IoMJvZEQVFez1svBg1yvq6OXe+hQ6nv0rVpx2kFfUS8rM6E3IoMJvZEQOHOnOsK2\ndy+9XCD1hb5rV+bUe3r0x49zFLAJvREJTOiNhMBX6Ovq3EKf6qGbbt1Y48ZT1HftYuzes4PaMJqK\nCb2REPgT+v79uZzKHn1WFoW+Qwdvj77UNdOyefRGJDChNxKC6mqGMBwPtqUIfffuLO/g69GXlXGe\nWRN6IxKY0BsJQXU10LdvywvdXHQR8NRTDT36sjL/ufWG0RRM6I24c/w4p9Tr08e7M7ZHD867msoe\nfVoap0rs2LFh6GbIEBN6IzKY0BtxZ/duCl1WFoXt2DHWoW/fnuGcVBZ6hw4dGoZuhg41oTcigwm9\nEXd27QK6dHHHqffto/ClpQHnnMMc81TH16N3hN6yboxIYEJvxIXaWuDVV7lcXe0t9J5T6s2dy0lH\nUp2MDIawDh/melmZhW6MyGFCb8SF5cuBu+/mcjChbymIuL//ihUs/TBwIEcGHz0ab+uMZMeE3ogL\ntbWccEPVLfSdOjFUsXdvyxN6gOGqv/wFmDEDmDYNyM1tmI1jGE3BhN6ICzU17HCtrDSP3qFjR+Cf\n/wQefJDvrVo1zK83jKZgQm/Ehdpavm/b5l/oO3WKq3lxoUMHToA+bZq7rVMnE3qj+ZjQG3Ghpobv\nW7cCGzcCAwaYR9+xIzB+PJCT491mQm80l5BCLyLPiEiFiKz1aJstIiUissr1Osdj2x0isklENorI\n1GgZbiQ3tbWcPm/bNmDlSgpcSxf6Dh2Ac8/1buvY0VIsjeaTHsY+fwfwZwDPebQpgEdV9VHPHUVk\nOIBLAAwH0BPAeyIyWFWPR8heI0WoqQHGjgU+/ZQDpgYM4NypdXUttzP25z931/dxMI/eiAQhPXpV\nXQyg1s8m8dM2A8ALqnpEVbcB2AxgQrMsNFKS2loK/fz5wMknc3CUI2olJcw4aWmcemrDMQONFfpD\nh4AFCyJrl5H8NCdG/yMRWSMifxORLFdbDwAlHvuUgJ69YXjhePQHDgDjxrHNCVO88w4wZUp87UsU\nAnXGDhniLu3syerVwFVXRd8uI7kIJ3TjjycA/Nq1fC+ARwBcG2Bf9dc4e/bs/y0XFBSgoKCgiaYY\nyUhtLTBmDJcdoW/fnh5pfT0wcmT8bEsk/Hn0hw8DX37J1NQuXby31dZyZqo9e9zz0frjvPOAe+7h\n3+C554BrrwWWLmXn+KWXRv57GE2jsLAQhYWFzT5Pk4ReVSudZRGZC+C/rtVSAL08ds13tTXAU+iN\nlkdNDdCzJ0d/nnIK20TYIXn22Vw2/At9eTnfa/0EVJ1spg0bgIkTA593+3Zg8WKWXbjuOnYCz53L\n403oEwdfJ3jOnDlNOk+TQjci4hlJvBCAk5EzD8AsEckQkX4ABgFY3iTLjJTl4EFWqGzXjqmVvXu7\nt3XsyEJmBvGXdbNzJ9937264vyP+GzYEP+/evewIX7qU6++8A7z7LmvsGKlHSI9eRF4AcDqAriJS\nDOAeAAUiMhoMy2wF8D0AUNUiEXkJQBGAowBuVFW/oRuj5VJbC3TuTK+9VSvvbXffbULviT+P3hHj\nQEKfkRGe0K9YwRIUBQXAH//ImL89SaUmIYVeVf09yD0TZP8HADzQHKOM1Kamhjn0/vj+92NrS6Lj\nKfRVVXwKCiX048YFF3pVd93/998HXnmFYn/55cBLL3G7CPD66/TyH3884l/LiDE2MtaIOY5Hb4Sm\nXz9g3TqK7xVXAE88QaFv1cp/jL62Fpg0KbjQ798PtGnDtFYAmDyZnbIzZrCPxMnm+eILvozkx4Te\niDnBPHrDmzFjKPJvvsn8+M8/Z4x+wIDgHn1JCftCHH74Q6C4mMt79zJtc9w4dtiK0HO/6CJO3+g8\nMezc6e74NZIbE3oj5phHHz4iwGWXMf3xxBMp9GVlwPDhgYU+J4cjbDdudLfPnw8UFXHZEfrrrgPu\nuINtXbty0Fpengl9KmJCb8QMVY7+XLTIPPrGcNlljM/ffz9DKTt2BBb6mhreREeP5uAph+pqoKKC\ny47QDx3Kv4cnPXq4s3p27uT5nFmvjOTFhN6IGdu2AZ98AjzzjHn0jWHQIIZuzj0X6NWL8ffhwwPH\n6Dt3Zshn1Sq2HTlCcfcVen/4evQAB2YZyY0JvREzPvuMIzLPOss7d94IzXnnsQN2xAggPZ3iHyh0\n4yv0zg0hHKH39egHDrTwTSpgQm/EjJUrmenx9ttWj6WpjBhBrzs7u6HQHznCDtiOHSn0a9Zw5KuT\nRROu0JeVMf1S1YQ+VTChN2KGI/RpaXwZjWfkSIpxVlbD0E1tLdtFWAMnM5O1axoj9E7oprycy3l5\nJvSpQFOLmhlG2Hz8MUVn5UpWrDSazvnnsxM1M5OFy5zBTYB7MheH0aMZvsnIYGzfEexgUzU6oZud\nOyny3bu7bxBG8mJ+lRF1HnqIg3jS0ykkRtPJyGBHbOvWQNu2wL597m2+aaujRzMds7qax4Qbutm9\nG1i71i305tEnPyb0RtTZsgW4+WbgmmuslkokycryjtM7qZUOAwYAX31FoR8yhNuPHQsu9BkZrDX0\n+OMU+txcE/pUwITeiCrHjzNOfOutwANWASmi+MbpfT36vn2Z0lpdzUFUWVnArl3BhR4AZs7k4Crz\n6FMHE3ojquzcySyQljgHbLTp3Nnt0d93H/tAPIW+Xz/eZGtq2DnrxNtDCf0553ASGIvRpw4m9EZU\nOHqUHuaWLQwhGJHHCd1UVlLof/97b6Hv2ZMefFkZhd4Jw4QS+nbtgDvvBMaPp9CXlXnXxPc3UMtI\nbEzojajw1lvMEPnqK9ZdMSKPI/QLFtALf+EFDqxyaNUKyM9n5k12NoU+HI8eoNAPH879rr6aov/a\na8C99/I8/uarNRIXS680okJZGcsdjBhhHn20yMlhOYSyMmDqVGDWrIb79O0LfPCB26MPV+gdRIA/\n/xn4z3/4rspc/uXLbYKYZMI8eiMqVFZSFJ5/3oQ+Wtx0E/D008B//wtMm+Z/n379+N6lC8tOfPll\n44TeYeZM3jAWLuRnOVMQGsmBCb0RFaqqOAy/vt6EPlr07QvcdRfQrVvg8FjfvnzPzgYuuAB49VUO\ntGqs0HsycaIJfbJhQm9EhcpK4Lvf5aO/CX30+OlPg4tuv37ACSewg7VvX9a0F2FbUznlFGDZMqbO\nGsmBCb0RFSorKSqrVzNzw4gOIsFLPvfty7CNM1Dtiiua580D7Bvo0sWmGUwmrDPWiApVVRSEESPi\nbUnLZswYZso4XHwx8+qby/jxzNsfNqz55zKij3n0RkRxZiOqrGTs2Igv7dqx9IRDp07AL37R/PMO\nGMBRt0ZyYEJvRIxjx1gUq66Og2q6do23RUa06N2bUxoayYEJvRExdu/mQJqPP6bnmG6BwZTFhD65\nMKE3IoYzWnLhQsbnjdTFhD65MKE3IobTybdwocXnU51evSj0qvG2xAgHE3ojYtTUMFyzcqV59KlO\np06c/MSzwNnGjZHp6DUijwm90SRU2fnqSU0NMG4cB9KYR5/6+IZv1q8H3n8/fvYYgQkp9CLyjIhU\niMhaj7ZsEVkgIl+KyLsikuWx7Q4R2SQiG0VkarQMN+LLc88BZ57p/eheU8M5Ydu1M4++JeAIvTNC\ntro6Mjn6RuQJx6P/O4CzfdpuB7BAVQcDeN+1DhEZDuASAMNdxzwuIvbUkIK8/TaH3j/3nLutupqe\n/PDh5tG3BHr3ZimE/v2ZUltdbbXqE5WQIqyqiwH4/vmmA3jWtfwsgAtcyzMAvKCqR1R1G4DNACZE\nxlQjUVAFFi0C/vEP4Lbb3B5dTQ2LZ33nOxw5aaQ2vXsDv/sdsH07JzSpqWHBNN+QnhF/mupt56qq\nM8FYBYBc13IPACUe+5UA6NnEzzASlC1bWDvl4ouBtm2BzZvZ7gj9zTez8JWR2vTuzb//0KEcCe2k\n13pOWG4kBs0e0qKqKiLBkqz8bps9e/b/lgsKClBQUNBcU4wYsWgRMHkyxX7sWOCzz4DBg91Cb7QM\npk8HBg4E7r+fE5o4Ql9by6JnRvMpLCxEYWFhs8/TVKGvEJHuqlouInkAKl3tpQB6eeyX72prgKfQ\nG8nFokXA6adz2RH6WbNM6Fsa7dszRJeb6+3RW4ds5PB1gufMmdOk8zQ1dDMPwFWu5asAvO7RPktE\nMkSkH4BBAJY38TOMBGXrVnrwAKsjrlrF5epqE/qWSE6OW+h79vTfIXv8uNWvjyfhpFe+AGAJgCEi\nUiwiVwN4EMBZIvIlgDNc61DVIgAvASgCMB/Ajao2di7V2LOHE1MDbo9elZ6cPbK3PDyFfuBA/x79\n7NnAY4/F3DTDRcjQjapeGmDTmQH2fwDAA80xykhs9uwBMjO53L07Zyvato1zkWZlBT3USEFychjO\nq61l+WJ/Hv2mTcChQ7G3zSBWX9BoNJ5CDwAnnwzMnw907Ai0ahU/u4z4kJtLIW/blsv+hL6srPkz\nWxlNxwYzGY1ClZ675z/tRRcBjz9u8fmWSk4OpxXs0oW/AX+hm7IyG0wVT0zojUaxfz9DNa1bu9u+\n9S3m1pvQt0xycoCDByn0nTt7C/qBA3QOSkstvz6emNAbjWL3bu+wDUDvfsYM64htqWRnA2lpfPf0\n6AsLgW98g6G+AwdM6OOJCb3RKHzj8w4//jEwbVrs7THiT1oaaxv5evTvvgt8/jlQUtKwpHEiMXs2\nQ4+pjHXGGo0ikNBPmsSX0TLJzXXH6B1BLyykJ//xx8CQIayHEy6LF3O07be/HRVzvVi/nr/rVMY8\neqNRBBJ6o2WTk+P26GtqgH376M2ffjrw3nusaLp7d/gzUr3yCvDnP0fXZofi4sbdhBrLxo3Af/4T\nvfOHgwm9EZJ//pNliVVN6A3/5OTQm3dCNx9/zLTbk0/mZCT9+jF8s39/eOdbt45lsA8ciK7dAIV+\n587onX/uXO9y3uFy003Az3/OjKXmYkJvhOSRR4Crr2Ys04Te8MeddwKXXML6N0eOAA8/zI7YE0+k\n8PfowZtAuB2y69YBeXnAJ59E1+7Dhyny0fTo33qraaGhf/wDWL0aeOqp5ttgQm+EpKIC+NWvgE8/\nNaE3/HPiiRwlLQJMmQKcdhpw660M2QAU+qys8Dpkq6o4inbWLMb5o0lZGdChQ/Q8+q1bGbppbMbR\noUO8Cc1gAwy7AAAXXklEQVScySeO5mJCbwTl+HH+440fzwkmTOiNUMyfD8yZQwEdNoxtPXu6PfpX\nXw0eklm3DjjpJD4RRFvoi4v5WYcOAfX1kT//W2+xpHdjPXqnQGCvXib0RgzYvZuP44MGuYXe6tkY\n4ZKZCUyYwOkGHY/+2muBf/878DGO0J9yirtgXrTYsYMTqHTvzifXSPPBB5ygx1Po33iDN8JgVFez\nc9uE3ogJFRVMncvMBNLT+ShqHr3RGJYt428oK4ulEnbvBv7+d/f2ffvcy3v3AmvXUuizsoCMDGDX\nLm5T5dNCJCkudgt9NMI327ezwmtdHZ+Oi4qAK68EXnwx+HG+Qt/cm50JvRGUigpmVABAnz78JzSh\nN5pC587sXJ0wgYK/aRPbTzsNWLKEseysLODpp4HRo7mtXz9WRgVYZuPcc4GPPgr9WRs3Ar/9LZf3\n7mUHsT+KiymmeXmR65AtK2N/FsDBYn368Km4rg647z7g9tvpMB0+HPgcjtB36sR+j+bm+ZvQG0Gp\nrKQ3BvAHW1xsQm80jawsCv3IkexofeklTiS+YQOwYAFH0l5zDfPwncF3ffu6hX79elbIvPfe0J/1\n6qvAL3/JUNH55wN//av//Ryh797dv9D/4Q+NH9H7zDMU9MOH+V1ycvjdd+/mTWDiRD5FOHMt+8MR\nepHIhG9M6I2gOKEbgEIPmNAbTaNzZ4rp0KH06teuZbGzw4cZy/7gA2bsdO5MgQPo0W/dyuWiIuD6\n6+mtn3ce8K9/Bf6sZcvoRd9yC0fZOrOg+eIZo/cN3Rw5wrTR1au9255/Hli4MPBnL11KES8r45NC\nq1b8n9mzh45Tt27MRioqch9z9Kj3ORyhB0zojShSWcn4om/oBjChN5qG04k/bBjFfsMGhmNOPhlY\nuZIZNmec4X2Mr0c/ZgwnORk3LrDQq1LoH3iAuegzZ3KUri9ffcXQyqBB/kM3a9YwO6ikhOd86ilO\noXn//Tx3oM9eupQ3px07gPx893d3hD4nx1voH36Y6acvvUQ7Dx5sKPQlJYGvaziY0Bt+mT0buO02\n79BN7958N6E3moKv0G/axFj9yJEU8Px892/NoW9fb49++HA6HOef39ADP34c+MUv6PEDwHXXAd/7\nHvD73/PYY8e897/nHuBHP+KEOf48emewVmkpj//VrzhKfNmywKN2N2/mk0RmJuP0jtBnZlK89+yh\ngDtC/957tO8vf+HAxClTgCeeMI/eiAFHjwIvv8wOMn+hG5spyGgKnTsDbdrQYWjfniGM99/nPLPn\nnMOXL05n7LFjFHAnLz8vr6Ewb9rEDthZs5ia2aoV8OSTzOHPyaEH77B2LfsFbrmF6z16NBTTJUt4\nAyopoYCPG8e+g8xMYMQIlnm48krvEMzSpYzBDxgAfPiht9Bv2cKbXatWFPolS4DLL+eTycyZvIHc\ndx+fXDznX87PN6E3osCHH9KLOHCAecxO6KZfP/5QPScdMYxw6dOHg6Cc6SaHDgXeeYdCf/vtwEMP\n+T9m+3Z69d260fsG6HxUVfEGMHcuc++XL2cRta1b2QfgyYgR3uGbu+/mE6tzvhEjKMR799KbLi2l\nEM+cyeWvvuJYAIczz+TI3+efZ9/CsWPAo49y3RH6RYt4kwH4f7Npk/t/acgQ3kB+9jNeE4ehQ3lD\nM4/eaDYPPcQfdCBefJFe0aRJ/DE6Hn1ODr0Nw2gK/ftzpKjDsGFMORwwgDXt0/yoUYcOfL31Fsss\nOLRuzSeEqioK/dy59IinT+eApKuv9j7PiBH04gGGZFatAn7wA/f2E06gx75kCXDzzbR1/36KcEmJ\nf6Ffs4azq61Zw5vIb3/LENA3v8nvtGePt0fvKfTt2jG0c+ut3nb6E/qRI+lwNWdydatH38I4fJix\nxsmTga99reH2fftYInb1asY8583zjpv26BE7W43UxgnDDBgQfL9+/fibffpp73YnfLNlC99zcoBL\nLwVOPbXhOUaOdHfePvggcNddDCN5cvrpfMKYN4/iXV/PczpCf9ZZ7n0nTQJef51ifMstvEGcdx5v\nOJ7fyVfoJ050n+PkkxvamZPDp4OtW91C3707b1QLFgS/TsEwj76FsW4dxb6khP8c3/mO9/bnnwcK\nCvi4OGkS85bbt4+LqUaKM3Qo0LVr6JIazz/POP3Mmd7teXnszK2vp4e/Zg1j6v446ywOtFq6lOmW\nl1/ecJ/Jk9kROm4cbRs7liJbXc3P8fTo09M5febIkXzKXbzYe+IdX6HPymL4xfHoAyHCzz50yHsO\n5m9/u3k17U3oWxgrV/K9pISPmy+84E4rUwUee4yZCAA7tP70J3dOs2FEklNOCTyQyZPBg/3fDPLy\n2CE6YABw4YUsm9Cunf9zZGdzMNb06aw94895mTiR3vSsWe629HSK85YtfLLwpVMnPvHOm+ct9AMH\nsi8iL4/rTqZaKKEHKPSdOnn3hX3rW8B//xv62ECY0LcQqqo4Mm/FCmY9lJS409Y++IDvCxdS1E8/\nnesZGUxRM4xo0KYNcNFFTT8+L4+e9IABTKO8++7g+99yC0OT11/vf3v79pw79uKLvdt79qSYB3qy\nHTWKT76DB7vbunblE4Yj1s6Nqlu30N9r6FB32MbThkB2h4PF6FsIP/6xe8DG9OnuCZt79WIu73e+\nw6nbbrrJPHgjOcjL41PpWWfRgx44MPT+ZWXBQ0X+xDQ/350p5I9Rozhi1rcz2bPzuLEeva/QA0yi\nePjh0Mf7w4S+BXD0KOuItG/PuPz99wO//jW3XXMNa3Ns3850sOefj6+thhEueXkMN4bqzPWkKSW2\ne/akxx6I664L/WTSGKGfOtUd8okUFrppAXzyCfORn3ySccTBg+nRb9sGTJtGD/7EE4Ebb2Qqm2Ek\nA44YNkbom8I3vgGcfXbg7fn57JQNRmNCN23bNhwH0FxEm1HoWES2AdgL4BiAI6o6QUSyAbwIoA+A\nbQAuVtXdPsdpcz7XCM2WLRTt3FwORklP56g7VWbddOzIDp9169jWpg3zkg0jWfjqK4r8li3eGTGJ\nyP79/H90Zo5qKiICVW10cLW5Hr0CKFDVMarq3INuB7BAVQcDeN+13uJ4+unwZ7yPBvfdx5g7wMEm\n553HZREODsnKon25ufSMTOSNZCMvj6EQpwZTItOuHQdoxWt2tuZ69FsBjFPVao+2jQBOV9UKEekO\noFBVh/ocFzePvr6ePe/hxMqaQ34+q9F5plzFkilTmDXzr38xbFNTQ6/eYexYljjYsCE+9hlGJDh6\n1Pt3nerE06N/T0RWiIjTX52rqs7sixUAcv0fGh+efZYV7qKJKqc/27Ejup8TjOJiplJ+8gkHgPj+\nM+TnszKgYSQzLUnkm0NzL9OpqrpTRLoBWODy5v+HqqqIJFQwvrSUKYbRpL6eI9tiKfT33MPyBN/7\nHm80xcXs1Pn3v/0/VTgj9gzDSH2aJfSqutP1XiUirwGYAKBCRLqrarmI5AHwK6uzZ88GwPzTKVMK\ncMYZBc0xJWx27nRPNuygyhGi7doxlh2sOqMjounpgeu+OOePxOzt4fLRRyzUdPHFfJxt25Zzcb74\nIvDaaw33v/BCy5c3jESnsLAQhYWFzT+RqjbpBaAdgI6u5fYAPgYwFcDDAG5ztd8O4EE/x6rDJZeo\nfutbqkePakw45xzVvn2923buVG3fXnXwYNUnnwx+/Jtvct+MDNUdO/zvs2KFKqA6fbp3++LFqu+8\nE56djzyiumSJanW16p13qh47Fnz/vn1VzzhD9Re/UF25UnXkSNU5c2hHdXV4n2kYRmLj0s5G63Vz\nYvS5ABaLyGoAywC8qarvAngQwFki8iWAM1zrAW4yHH6/bRur0zWVvXs5iW84+PPo169nJblbb2UV\numC8+y6HWp9xRuB5KJ0UKt/QzcsvM5TiyZEjHLDka8999wEXXMD83YcfZrw9EEePcsTfI4+w8FFx\nMTMRxo/nBAfNSecyDCP5abLQq+pWVR3tep2kqr9xtdeo6pmqOlhVp6pPDr0nX3zBEMNjj4VXgvNn\nP6Og+bJgAXDHHQ2nCvNHeTmzbg4edLc5U5RNmhRa6D/4gCI/ciRrWfhj1y5W0fMN3WzZ0lD8V60C\nrr3Wu/2uu/j6/e9p0y23sGiSJ/fdB8yfz+XiYqZJjhzJejbLl7O0wbRpLLtqGEbLJq4jYz/6CPj6\n18Ob/La2ljO4/O53DbctXEjh3rzZ/7EVFaySd+wYRTgnh163w/r1HBk6bBi3B+qsraignWPHsr5F\nIKGvrma9iv37vXPpv/qqodAvX853J45eUsLr8sMfsv7MH/5Az/6NN+ixjx/PzKHHHgO++13e+LZu\n5YCRtDTWmH/xRV7TtDTrdDUMI85Cv3gxhT43lwJ75EjgfZcupRj/4x/eIg1w9vjevf3P9P7llxxO\nfNNNFMTsbNaYrqpy7+MIfVoaS6cuWkSx9E31X7iQNavT0yn0/j4P4Hfp2tV7CrDjx/n5xcXu9EtV\nzjIzYwYn+wBYdvW007wnRZgwgTcfZ/Lgq68G/u//mGFzyy28gTglVE89lU8OyTCIxDCM2BA3oT9y\nhILqCGdODuPnDz3EWLYvS5bQs73oInqzDpWV9IIvu8y/8L7yCqf2GjWKYYy8PIqwE6dXpdAPH871\nSZPoKc+a1fAJ4bXX3LPMOHM++hv9Wl3tFnrHgy8vZ9mBDh1o8xlnsIDY8uXAnXfS9ooKfk/fdMhW\nrRiqefllzo5TWsqpzH76U456Xb/eLfTOsb16Bbz0hmG0MOI23OChh5jOONQ1ZjY/n8L57ruswfLt\nb3vvv2QJY/QDB9Ljvflmerhbt3J97Fh6uatXM4zjTNlVXs56GAcOAG+/TaHPzKTQb9nCjsy0NPdI\n2YsvprB+/jlDKIMGsX3VKt6Y/vY3rqenM9Tz4YcceOTcKACeu0sXetVvv83ywN27045Dhzgn5Pr1\nvAbFxbT9ggv4tPLxx8Af/9jwenmWT3WKOWVn09t/7jn3MRMm0DYTesMwHOLm0d9zDzBnjjuX2xH6\njRuB999n4S2A2Tg/+hFDHBMnsvJiQQGXFy5k6d0bbuCciqtWcYqwe+91f05FBUNDo0ezI7V7d3rb\nVVXAFVdw1OiJJ7rtGDqUHaGTJ1N0HW6/HfjlL72rO44ZwyeMggLviXud0M2oUfS4r7uOnnv//ixH\n8Oab7Dg9dIh2paezXvyf/sSSBP7mkgzEhReyvIHj0bdr5775GIZhAHEU+lGjvGs45+fTy929m2K+\nZAlF/4knWCu9Tx93muBdd9GrnTeP26dPp6dfWcmCXR995I73V1RQ3EeP5ohVJ3RTUUGv/emngZ/8\npKF9p53G8wAMzyxa1HC2pd/9juGmE0/0zopxZnD/6U/5Hc49lxk0AwbQy3/9dd5g5sxxP7mMHcvv\nMGpUw0mLgzFjBt89pzmbNMkGQxmG4SZuoZslS7xnZMnPpwAOGUJh/OtfKaK33cYOx7o6976jRtHr\n96RVK8bpv/99CvKnn1Lwysvp0TvZJ3l5FMFXX2W75/yQnpx0Eo+tqmJ4aOhQFgnzxKlEd/31vGE4\nkxc7Hr3DZZcxf37AAG7bvJmhp8su8z7fQw81vmxCfj5vQpGeqMAwjNQhbh69r9eanw8sW0ZBnTmT\nYZxJkxi2SUtzz9ASjLlz6SlPmeKeB9Xx6DMz6fXm5bH4/8cfB54xHuCN42tfo1e/bh2FPxAXXcS+\ngU8+4brj0TtMm0bhHzjQnQ0zdmzD80yc2HC+ynD4+tfNgzcMIzAJM8NUfj47RocOZbx98WLggQcY\nimksU6bQ4z90iIOjnFrrd91F8e7alX0AwYQe4MwyH3xAoR8xIvB+bdoAf/87wzBffMFMHs/Z6Fu3\n5k1s4kQKfVpa6BlpDMMwIkVCCT3gzsJpDqeeytBNeTm9dydEdO21nP/RCav486o9Oessjrpduza4\n0AMshnbrrcDpp/P8vh62M6Bp+HB27HreCAzDMKJJwlRzdipBRkLoMzMptkuWMGzjizNvYyiPftQo\nZrSUlAQP3TjcfDO9+3ffDbxPhw6cnNswDCNWJIxHn5HBgUNDhkTmfCNG0BvP9TPtSU4OP8vfTcCT\ntDSGgdLTwy8l8IMf+C8LbBiGES8SRugBerpNicn7Y8QI4L33/It5enr4XvXUqTyXdXYahpGsJEzo\nJtKMGOGu6tgcLruM5QoMwzCSlYTy6COJk9USKjwTiowMDtYyDMNIVlJW6AcPZlpjcz16wzCMZCdl\nhb51a4ZvrB67YRgtHVHfouux+FARjcXnOlUkrSPVMIxUQESgqo1WtJQWesMwjFSiqUKfsqEbwzAM\ng5jQG4ZhpDgm9IZhGCmOCb1hGEaKY0JvGIaR4pjQG4ZhpDgm9IZhGCmOCb1hGEaKY0JvGIaR4kRF\n6EXkbBHZKCKbROS2aHyGYRiGER4RF3oRaQXgMQBnAxgO4FIRGRbpz4kFhYWF8TYhLJLBzmSwETA7\nI43ZmRhEw6OfAGCzqm5T1SMA/g1gRhQ+J+okyx8/GexMBhsBszPSmJ2JQTSEvieAYo/1ElebYRiG\nEQeiIfRWltIwDCOBiHiZYhGZCGC2qp7tWr8DwHFVfchjH7sZGIZhNIGEqEcvIukAvgAwBUAZgOUA\nLlXVDRH9IMMwDCMs0iN9QlU9KiI3AXgHQCsAfzORNwzDiB9xmWHKMAzDiB1RGxkrIm1EZJmIrBaR\nIhH5TYD9/uQaWLVGRMZEy56m2igiBSKyR0RWuV53x9JGH1tauWz4b4DtcbuWPnYEtDNRrqeIbBOR\nz102LA+wT9yvZyg7E+h6ZonIyyKywfW/NNHPPolwPYPaGe/rKSJDPD57lcuWH/vZr3HXUlWj9gLQ\nzvWeDmApgNN8tp8L4C3X8ikAlkbTnibaWABgXqztCmDrLQD+6c+eRLiWYdqZENcTwFYA2UG2J8T1\nDMPORLmezwK4xrWcDiAzQa9nKDsT4nq6bEkDsBNAr+Zey6jWulHVetdiBhivr/HZZTp44aGqywBk\niUhuNG3yJQwbAaDRvdyRRkTywT/wXPi3J+7XEgjLTgRpjzXB7EiI6+ki1PWK6/UUkUwAX1fVZwD2\n06nqHp/d4n49w7QTSJzf55kAtqhqsU97o69lVIVeRNJEZDWACgALVbXIZxd/g6vyo2mTL2HYqAAm\nuR6R3hKR4bG0z4PfA/g5gOMBtsf9WroIZWeiXE8F8J6IrBCR6/1sT5TrGcrORLie/QBUicjfReQz\nEXlaRNr57JMI1zMcOxPhejrMAvAvP+2NvpbR9uiPq+polxGTRaTAz26+d8+Y9g6HYeNn4KPTKAB/\nBvB6LO0DABE5H0Clqq5CcG8jrtcyTDvjfj1dnKqqYwCcA+CHIvJ1P/vE9Xq6CGVnIlzPdABjATyu\nqmMB7Adwu5/94n09w7EzEa4nRCQDwDcB/CfQLj7rQa9lTMoUux6P/h+AcT6bSgH08ljPd7XFnEA2\nqmqdE95R1fkAWotIdozNmwRguohsBfACgDNE5DmffRLhWoa0M0GuJ1R1p+u9CsBrYI0mTxLheoa0\nM0GuZwmAElX91LX+MiioniTC9QxpZ4JcT4A39pWuv7svjb6W0cy66SoiWa7ltgDOArDKZ7d5AK50\n7TMRwG5VrYiWTU2xUURyRURcyxPAlFR/cfyooap3qmovVe0HPs59oKpX+uwW12sZrp2JcD1FpJ2I\ndHQttwcwFcBan93ifj3DsTMRrqeqlgMoFpHBrqYzAaz32S3u1zMcOxPherq4FHSW/NHoaxnxAVMe\n5AF4VkTSwBvK86r6voh8DwBU9a+q+paInCsim8HHqKujaE+TbATwbQA/EJGjAOpBAYs3CgAJdi39\n0cBOJMb1zAXwmuv/OR3AP1X13QS8niHtRGJcTwD4EYB/ukIOWwBck4DXM6SdSIDr6bqpnwngeo+2\nZl1LGzBlGIaR4thUgoZhGCmOCb1hGEaKY0JvGIaR4pjQG4ZhpDgm9IZhGCmOCb1hGEaKY0JvGIaR\n4pjQG4ZhpDj/H0alzV+vYglmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x104aa4208>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.errorbar(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# the normalised residual for the data\n",
    "def residual(p, just_generative=False):\n",
    "    v = p.valuesdict()\n",
    "    generative = v['a'] + v['b'] * x\n",
    "    M = 0\n",
    "    while 'a_max%d' % M in v:\n",
    "        generative += gauss(x, v['a_max%d'%M], v['loc%d'%M], v['sd%d'%M])\n",
    "        M += 1\n",
    "\n",
    "    if just_generative:\n",
    "        return generative\n",
    "    return (generative - y) / dy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create a Parameter set for the initial guesses\n",
    "def initial_peak_params(M):\n",
    "    p = lmfit.Parameters()\n",
    "    \n",
    "    # a and b give a linear background\n",
    "    a = np.mean(y)\n",
    "    b = 1\n",
    "    \n",
    "    # a_max, loc and sd are the amplitude, location and SD of each gaussian component\n",
    "    a_max = np.max(y)\n",
    "    loc = np.mean(x)\n",
    "    sd = (np.max(x) - np.min(x)) * 0.5\n",
    "\n",
    "    p.add_many(('a', np.mean(y), True, 0, 10), ('b', b, True, 1, 15))\n",
    "\n",
    "    for i in range(M):\n",
    "        p.add_many(('a_max%d'%i, 0.5 * a_max, True, 10, a_max),\n",
    "                   ('loc%d'%i, loc, True, np.min(x), np.max(x)),\n",
    "                   ('sd%d'%i, sd, True, 0.1, np.max(x) - np.min(x)))\n",
    "    return p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving with `minimize` gives the Maximum Likelihood solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[Variables]]\n",
      "    a:        1.51236732 (init= 10)\n",
      "    b:        10.5608914 (init= 1)\n",
      "    a_max0:   193.129348 (init= 141.2683)\n",
      "    loc0:     5.03878685 (init= 5)\n",
      "    sd0:      0.57126003 (init= 2)\n",
      "[[Correlations]] (unreported correlations are <  0.500)\n"
     ]
    }
   ],
   "source": [
    "p1 = initial_peak_params(1)\n",
    "mi1 = lmfit.minimize(residual, p1, method='differential_evolution')\n",
    "\n",
    "lmfit.printfuncs.report_fit(mi1.params, min_correl=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From inspection of the data above we can tell that there is going to be more than 1 Gaussian component, but how many are there? A Bayesian approach can be used for this model selection problem. We can do this with `lmfit.emcee`, which uses the `emcee` package to do a Markov Chain Monte Carlo sampling of the posterior probability distribution. `lmfit.emcee` requires a function that returns the log-posterior probability. The log-posterior probability is a sum of the log-prior probability and log-likelihood functions. \n",
    "\n",
    "The log-prior probability encodes information about what you already believe about the system. `lmfit.emcee` assumes that this log-prior probability is zero if all the parameters are within their bounds and `-np.inf` if any of the parameters are outside their bounds. As such it's a uniform prior. \n",
    "\n",
    "The log-likelihood function is given below.\n",
    "\n",
    "To use non-uniform priors then should include these terms in `lnprob`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This is the log-likelihood probability for the sampling.\n",
    "def lnprob(p):\n",
    "    resid = residual(p, just_generative=True)\n",
    "    return -0.5 * np.sum(((resid - y) / dy)**2 + np.log(2 * np.pi * dy**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To start with we have to create the minimizers and *burn* them in. We create 4 different minimizers representing 0, 1, 2 or 3 Gaussian contributions. To do the model selection we have to integrate the over the log-posterior distribution to see which has the higher probability. This is done using the `thermodynamic_integration_log_evidence` method of the `sampler` attribute contained in the `lmfit.Minimizer` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 310 (-4394.1394141818528, 0.78699415494884306)\n",
      "1 310 (-1038.5082602579189, 3.5440488907950112)\n",
      "2 310 (-923.39096623751232, 5.0783786979153547)\n",
      "3 310 (-923.6894069769952, 5.6418844327914712)\n"
     ]
    }
   ],
   "source": [
    "# Work out the log-evidence for different numbers of peaks\n",
    "total_steps = 310\n",
    "burn = 300\n",
    "thin = 10\n",
    "ntemps = 15\n",
    "workers=4\n",
    "log_evidence = []\n",
    "mini = []\n",
    "res = []\n",
    "\n",
    "# set up the Minimizers\n",
    "for i in range(4):\n",
    "    p0 = initial_peak_params(i)\n",
    "    mini.append(lmfit.Minimizer(lnprob, p0))\n",
    "    out = mini[i].minimize(method='differential_evolution')\n",
    "    res.append(out)\n",
    "\n",
    "# burn in the samplers\n",
    "for i in range(4):\n",
    "    # do the sampling\n",
    "    out = mini[i].emcee(steps=total_steps, ntemps=ntemps, workers=workers, reuse_sampler=False,\n",
    "                       params=res[i].params)\n",
    "    # get the evidence\n",
    "    print(i, total_steps, mini[i].sampler.thermodynamic_integration_log_evidence())\n",
    "    log_evidence.append(mini[i].sampler.thermodynamic_integration_log_evidence()[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we've burned in the samplers we have to do a collection run. We thin out the MCMC chain to reduce autocorrelation between successive samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 410 (-4394.2785416442975, 0.83244601315891487)\n",
      "1 410 (-1038.6450874311915, 3.5634191857072892)\n",
      "2 410 (-923.45302756653928, 5.0966495848754221)\n",
      "3 410 (-923.82736691825357, 5.6538341407954249)\n",
      "0 510 (-4394.3487587435875, 0.85588467626166675)\n",
      "1 510 (-1038.6934181814577, 3.5598495634503706)\n",
      "2 510 (-923.53824600518396, 5.1203799324430292)\n",
      "3 510 (-923.82578592837626, 5.6281725406512351)\n",
      "0 610 (-4394.401228206887, 0.87176207214633905)\n",
      "1 610 (-1038.7503716831113, 3.5601188237437782)\n",
      "2 610 (-923.61139504744551, 5.1481129345677346)\n",
      "3 610 (-923.77987084071901, 5.5961124574821497)\n",
      "0 710 (-4394.4249700306227, 0.87627960796453408)\n",
      "1 710 (-1038.7884107190625, 3.5548852825279482)\n",
      "2 710 (-923.77476076212815, 5.2038245396746561)\n",
      "3 710 (-923.83444327268512, 5.6099999289124298)\n",
      "0 810 (-4394.4437133300844, 0.88390942394744343)\n",
      "1 810 (-1038.8249727101168, 3.5588724286415072)\n",
      "2 810 (-923.92607632064164, 5.248356094005203)\n",
      "3 810 (-923.8713318333505, 5.6097200131257523)\n",
      "0 910 (-4394.4482814566245, 0.88343021249966114)\n",
      "1 910 (-1038.8457492992793, 3.5558249795581105)\n",
      "2 910 (-924.07527838691749, 5.3002800701458455)\n",
      "3 910 (-923.98092996478613, 5.625742358755133)\n"
     ]
    }
   ],
   "source": [
    "for j in range(6):\n",
    "    total_steps += 100\n",
    "    for i in range(4):\n",
    "        # do the sampling\n",
    "        res = mini[i].emcee(burn=burn, steps=100, thin=thin, ntemps=ntemps,\n",
    "                         workers=workers, reuse_sampler=True)\n",
    "        # get the evidence\n",
    "        print(i, total_steps, mini[i].sampler.thermodynamic_integration_log_evidence())\n",
    "        log_evidence.append(mini[i].sampler.thermodynamic_integration_log_evidence()[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x1096be7b8>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEPCAYAAACOU4kjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucVXW9//HXO1EEJdS0QMTUnBTUEMnoojjnx09/2vkJ\nmBaYmReyAhXtYCni0SnzQl5Qj6nV0QRKDna8YRJekCntiJSAkugJ85IgdhIFzBu3z/ljrYHtOMze\nM7P3XnvveT8fj/2Ytb/ru/f6LJfMZ77rs9Z3KSIwMzMrpg9lHYCZmdUeJxczMys6JxczMys6Jxcz\nMys6JxczMys6JxczMyu6iksukhokLZO0MH0dlbNugqSlkp6VdERO+yBJi9N112YTuZmZNam45AIE\ncHVEDExfvwGQ1B8YCfQHjgRukKT0MzcCoyOiDqiTdGQWgZuZWaISkwuAWmgbDkyPiHUR8SLwHDBY\nUm+gR0TMT/tNBUaUJ0wzM2tJpSaXMyU9KelmSTukbbsCy3L6LAP6tNC+PG03M7OMZJJcJD2Y1kia\nv4aRnOLaEzgQWAFclUWMZmbWfl2y2GhEHF5IP0n/Dtybvl0O9M1ZvRvJiGV5upzbvryF7/IkamZm\n7RARLZUqWlVxp8XSGkqTY4DF6fJMYJSkbSTtCdQB8yPiVWCNpMFpgf9E4O6WvjsiavZ10UUXZR6D\n98/719n2rTPsX3tlMnLJY5KkA0muGnsB+BZARCyRdDuwBFgPjI3Nez4WuBXoBsyKiNllj9rMzDap\nuOQSEV9vZd2lwKUttD8BHFDKuMzMrHAVd1rM2qe+vj7rEErK+1e9annfoPb3r73UkXNq1URSdJZ9\nNTMrFklELRT0zcys+jm5mJlZ0Tm5mJlZ0Tm5mJlZ0Tm5mJlZ0Tm5mJlZ0Tm5mJlZ0Tm5mJlZ0Tm5\nmJlZ0Tm5mJlZ0Tm5mJlZ0Tm5mJlZ0Tm5mJlZ0Tm5mJlZ0Tm5mJlZ0Tm5mJlZ0VXcY47NrLqtWwdr\n1sDq1cnrH/8AP6ev88kkuUj6MtAA7AscHBELctZNAE4FNgDjIuKBtH0QcCuwLTArIs5K27sCU4GD\ngJXAyIh4qWw7Y1YjIuDtt9+fGFavfv/7QpbXroUPfxh69kx+br89bLVV1ntn5ZbVyGUxcAzwk9xG\nSf2BkUB/oA/wkKS69PnENwKjI2K+pFmSjoyI2cBoYGVE1EkaCUwCRpVzZ8yytmHD5l/wW/qlny85\nrFkDW2+9OTE0JYfmy337wv77b7lf9+6gNj8U1ypVe49lJsklIp6F5NnMzQwHpkfEOuBFSc8BgyW9\nBPSIiPlpv6nACGA2MAy4KG2/A7i+xOGbFU0EvPde20cKzd+/8w706PH+X/ItJYc+fbacND78Ydhm\nm6z/i1itqLSay67AvJz3y0hGMOvS5SbL03bSny8DRMR6Sasl7RQRr5chXuvENm6EN99s/0ihaVlq\nfaTQsyf06gX77LPlfttvDx/y5TlWQUqWXCQ9CPRqYdX5EXFvqbbbmoaGhk3L9fX11NfXZxGGVYC1\nazs2UlizJilUd+/e8i/83Pd1dVtOGj17QteuWf/XMNussbGRxsbGDn+PIsPLOCTNBcY3FfQlnQcQ\nEZen72eTnPJ6CZgbEf3S9uOBIRExJu3TEBHzJHUBVkTELi1sK7LcVyuOCHjrrY4VnFevTmoUrY0U\nWksGTe979HCh2mqfJCKizZWXSjgtlhv0TOA2SVeTnO6qA+ZHREhaI2kwMB84Ebgu5zMnkZxOOw6Y\nU7bIrSzGjYOZM5Ok8OabyV/6+RLAHnu0nii6dXPR2ayUsroU+RiS5LAzcJ+khRFxVEQskXQ7sARY\nD4zNGW6MJbkUuRvJpciz0/abgWmSlpJciuwrxWrIQw/Br3+d/Nxpp2S0sPXWWUdlZvlkelqsnHxa\nrPqsXQuf+hT86EcwbFjW0Zh1Tu09LebrS6xiXXst7LUXHH101pGYWVt55GIVaflyGDAA5s2DvffO\nOhqzzqu9IxcnF6tIX/1qMmr54Q+zjsSsc3NyycPJpXo0NsJJJ8GSJbDddllHY9a5ueZiNWHdOjjj\nDLj6aicWs2rm5GIV5cc/ht694UtfyjoSM+sInxazivHqq8lsu48+Cvvum3U0ZgauueTl5FL5TjoJ\nPvax5L4WM6sM1Tz9ixm//z3MmQPPPJN1JGZWDK65WOY2bIDTT4crr0ymdzGz6ufkYpm76SbYYQcY\nOTLrSMysWFxzsUz9/e/Qvz/MnZsU882ssrign4eTS2X6xjeSU2GTJ2cdiZm1xAV9qzqPPw6zZrmI\nb1aLXHOxTDQV8SdNSh7eZWa1xcnFMnHzzbDttvC1r2UdiZmVgmsuVnYrVyZF/PvvhwMPzDoaM2uN\nC/p5OLlUjjFjYKut4Prrs47EzPJxQd+qwhNPwF13uYhvVusyqblI+rKkpyVtkHRQTvsekt6RtDB9\n3ZCzbpCkxZKWSro2p72rpBlp+zxJHy/3/lhhNm5MptO/5BLYccesozGzUsqqoL8YOAb4XQvrnouI\ngelrbE77jcDoiKgD6iQdmbaPBlam7ZOBSaUM3NpvyhSIgFNOyToSMyu1TJJLRDwbEX8utL+k3kCP\niJifNk0FRqTLw4Ap6fIdwNCiBWpFs2oVTJiQ1Fk+5GsUzWpeJf4z3zM9JdYo6ZC0rQ+wLKfP8rSt\nad3LABGxHlgtaaeyRWsFufBCGD4cPv3prCMxs3IoWUFf0oNArxZWnR8R927hY68AfSPijbQWc7ek\n/YoVU0NDw6bl+vp66uvri/XV1oonn4T/+A8X8c2qQWNjI42NjR3+nkwvRZY0FxgfEQtaWw+sAB6O\niH5p+/HAkIgYI2k20BAR8yR1AVZExC4tfJcvRc5ABAwZAiecAN/+dtbRmFlbtfdS5Eo4LbYpaEk7\nS9oqXd4LqAOej4gVwBpJgyUJOBG4J/3YTOCkdPk4YE7ZIre8fvlLePttOO20rCMxs3LKZOQi6Rjg\nOmBnYDWwMCKOknQs8H1gHbARuDAi7ks/Mwi4FegGzIqIcWl7V2AaMBBYCYyKiBdb2KZHLmW2Zg30\n6wd33AGf/WzW0ZhZe/gO/TycXMpv/Hh44w245ZasIzGz9vId+lZRnn4apk5NfppZ51MJNRerMRFw\n5pnJ5ccf/WjW0ZhZFpxcrOhuvz2Z+XjMmKwjMbOsuOZiRfWPfyRF/Ntug0MPzToaM+soF/TzcHIp\nj/POg+XLYdq0rCMxs2JwcsnDyaX0/vu/4QtfgMWLoXfvrKMxs2Ko5psorQZEwLhxcP75Tixm5uRi\nRXLXXbBsWXKVmJmZT4tZh739NvTvDz//OfzTP2UdjZkVk0+LWWYuuyyZ3sWJxcyaeORiHfLcc0li\nWbQIdtst62jMrNg8crFMnH02fPe7Tixm9n6eW8za7de/hqVL4c47s47EzCqNk4u1y7vvwllnwY03\nwjbbZB2NmVUanxazdvnRj+DAA+GII7KOxMwqkQv61mYvvgiDBsGCBfDxj2cdjZmVkgv6Vjbf+U7y\ncmIxsy1xzcXaZPbsZO6w6dOzjsTMKplHLlaw995L5g+79lrYdtusozGzSlZQcpHUXdI+xdqopCsk\nPSPpSUl3SuqZs26CpKWSnpV0RE77IEmL03XX5rR3lTQjbZ8nySdrSmTyZNh3X/jnf846EjOrdHmT\ni6RhwELg/vT9QEkzO7jdB4D9ImIA8GdgQvrd/YGRQH/gSOAGSU2FpBuB0RFRB9RJOjJtHw2sTNsn\nA5M6GJu14OWX4cor4Zprso7EzKpBISOXBmAw8AZARCwE9urIRiPiwYjYmL59HGi6v3s4MD0i1kXE\ni8BzwGBJvYEeETE/7TcVGJEuDwOmpMt3AEM7Epu1bPx4OP102KtDR97MOotCksu6iFjVrG1jiz3b\n51RgVrq8K7AsZ90yoE8L7cvTdtKfLwNExHpgtaSdihhfpzdnDvzhD3DuuVlHYmbVopCrxZ6WdALQ\nRVIdMA74r3wfkvQg0KuFVedHxL1pn4nA2oi4rQ0xt1tDQ8Om5fr6eurr68ux2aq2dm3yjJbJk6F7\n96yjMbNSa2xspLGxscPfk/cmSknbAROBpuL6/cDFEfFuhzYsnQycBgxt+i5J5wFExOXp+9nARcBL\nwNyI6Je2Hw8MiYgxaZ+GiJgnqQuwIiJ2aWF7vomyHa66Ch56CGbNArX5Niozq3btvYkykzv002L8\nVcBhEfFaTnt/4DbgMySnux4C9o6IkPQ4yahpPnAfcF1EzJY0FjggTTSjgBERMaqFbTq5tNErr8Cn\nPgX/9V/wyU9mHY2ZZaFkyUXSQ8BxTXWXtJ4xPSL+X7siTb5jKbAN8Hra9FhEjE3XnU9Sh1kPnBUR\nTVepDQJuBboBsyJiXNreFZgGDARWAqPSiwGab9PJpY1OOCG5C//SS7OOxMyyUsrksigiDszXVumc\nXNrmd7+Dr30NnnkGttsu62jMLCulnFtsQ+6NiZL2oLhXi1mFWb8ezjgjqbc4sZhZexRytdhE4BFJ\nv0vfDwG+WbqQLGs33AAf/Sgcd1zWkZhZtSqooC9pF+CzQADzcovw1cKnxQrzt7/B/vsnp8X69cs6\nGjPLWkmvFpPUB9iDZKQTABHxu9Y+U2mcXApzyimw885wxRVZR2JmlaC9ySXvaTFJk0jm+1oCbMhZ\nVVXJxfJ77DF44AF49tmsIzGzaldIzeUYYJ+IeK/UwVh2NmxI5g674gro0SPraMys2hVytdhfSO5J\nsRr2058mSeX447OOxMxqQSEjl3eARZLmAE2jl2i6idGq32uvwUUXJRNUeooXMyuGQm6iPDldbOoo\nkuQypeVPVCYX9Lfsm99MJqX0s1rMrLlSXy3WHdg9Iqq21Ovk0rI//AGGDUvuxN9hh6yjMbNKU7I7\n9HOeRDk7fV+MJ1FaBdi4MSniX365E4uZFVcmT6K0ynDLLdClC5x4YtaRmFmtKaSgvy4iVun9lV7P\nLVblXn8dJk6E2bPhQ4X8iWFm1gYlexKlVbZ//ddk7rCBA7OOxMxqUWZPoiw3F/Q3W7gQjjoKliyB\nnXbKOhozq2RV9STKLDi5JDZuhEMOgVNPhW98I+tozKzSFX1uMUn35rwNkvtbNr2PiGFt3Zhlb9q0\n5Hktp56adSRmVstaq7lclf48BugF/IIkwRwP/K3EcVkJrFoF550HM2e6iG9mpVVIzeWJiBiUr63S\n+bQYnH02vPUW/OxnWUdiZtWilI857i7pEzkb2gvo3tYN5ZJ0haRnJD0p6U5JPdP2PSS9I2lh+roh\n5zODJC2WtFTStTntXSXNSNvn5T6S2TZbvBhuuw0uuyzrSMysMygkuXwHmCvpt5J+C8wFzu7gdh8A\n9ouIAcCfgQk5656LiIHpa2xO+43A6IioA+okHZm2jwZWpu2TgUkdjK3mRMAZZ8D3v588CMzMrNTy\n3ucSEbMlfRLYl6Sw/2xHn+0SEQ/mvH0cOLa1/pJ6Az0iYn7aNBUYQTIlzTDgorT9DuD6jsRWi6ZP\nhzffTCaoNDMrh9auFhsaEXMkHcv7rxb7RHoO7s4ixXAqMD3n/Z6SFgKrgQsi4lGgD7Asp8/ytI30\n58sAEbFe0mpJO0XE60WKr6q9+SZ897vwq1/BVltlHY2ZdRatjVyGAHOAo9k83X6uVpOLpAdJrjJr\n7vyIuDftMxFYGxG3peteAfpGxBuSDgLulrRfnn0oWENDw6bl+vp66uvri/XVFesHP4DDD4fPfz7r\nSMysGjQ2NtLY2Njh7ynkarEuEbG+w1v64PeeDJwGDN3S3f6S5gLjgRXAwxHRL20/HhgSEWMkzQYa\nImKepC7AiojYpYXv6nRXiy1ZAocdBn/6E3zsY1lHY2bVqJRXiz0v6aeShkrFeU5hWoz/LjA8N7FI\n2lnSVunyXkAd8HxErADWSBqcxnAicE/6sZnASenycSSjrU4vAsaNgwsucGIxs/IrdG6x/w+MAg4C\n7gVmRMQj7d6otBTYBmiqizwWEWPT+s73gXUkMy9fGBH3pZ8ZBNwKdANmNT1mWVJXYBowEFgJjIqI\nF1vYZqcaufzqV3DxxbBgQTKtvplZe5RlbjFJOwLXAV+NiKoqD3em5PLWW9CvH/ziFzBkSNbRmFk1\nK+VpMSTVS7oRWAB0Bb7S1g1Z+VxyCRx6qBOLmWWnkNNiLwKLgBnAvRHxjzLEVXSdZeTy5z8nV4Y9\n9RTsumvW0ZhZtSvZaTFJPSNidbsjqxCdIblEwBe/CEOHwjnnZB2NmdWCUp4W6yVpjqSn0w0NkHRB\nmyO0krvnHnjpJTjrrKwjMbPOrpDk8jPgfGBt+v4pkmn3rYK88w585zvwb/8GW2+ddTRm1tkVNCty\nRDze9CY9t7SudCFZe1x+ORx8cHJKzMwsa4XcAfF3SXs3vZF0HMkd81Yhnn8err8eFi3KOhIzs0Qh\nBf1PAD8FPgesAl4ATmjpRsVKVssF/WHD4HOfgwkT8vc1M2uLkt9EKWn7tP+bbd1IJajV5HLffUmt\nZfFi6No162jMrNaU9CZKgPT+lul5O1rZvPtucmXYddc5sZhZZSk4uaT65O9i5XLllXDAAXDkkfn7\nmpmVU1unNHTJuEK89BJMngxPPJF1JGZmH9SmiSurWa3VXI49FgYMgAsvzDoSM6tl7a255B25SFrM\n+x9zDMkjiP8A/DAiVrZ1o9YxDzyQXHb8y19mHYmZWcsKOS02G1gP3EaSYEYB3YG/kTxf5ehSBWcf\ntHZt8hCwa66BbbfNOhozs5YVklz+b0QMzHn/lKSFETEwHdVYGU2eDHvvDUc7pZtZBSskuWwlaXDT\nFDCSPsPmq8zWlywy+4Bly+CKK+Dxx/P3NTPLUiHJZTTw8/QmSoA3gdHp448vK1lk9gHnnANjxsAn\nPpF1JGZmrWvLHfo9Aar12S7VfrXY3LlwyimwZAl07551NGbWWZTsDn1JO0iaDDwMPCzpqqZE016S\nLpb0pKRF6bNi+uasmyBpqaRnJR2R0z5I0uJ03bU57V0lzUjb50n6eEdiq0Tr1sEZZ8DVVzuxmFl1\nKOQO/VuANcCXga+QnBb7eQe3+6OIGBARBwJ3AxcBSOoPjAT6A0cCN0hqypg3AqMjog6ok9R0X/po\nYGXaPhmY1MHYKs7118Nuu8Exx2QdiZlZYQqpuXwiIr6U875B0pMd2WizyS+3B15Ll4cD0yNiHfCi\npOeAwZJeAnpExPy031RgBMll0sNIkxNwB3B9R2KrNCtWwKWXwqOPgto8MDUzy0YhI5d3JB3a9EbS\nIcDbHd2wpEsk/RU4mc0XBuwKLMvptoxkPrPm7cvZPM9ZH+BlgIhYD6yWtFNH46sU3/sejB4N++yT\ndSRmZoUrZOTybWBqTp3lDeCkfB+S9CDQq4VV50fEvRExEZgo6TzgGuCUAmNut4aGhk3L9fX11NfX\nl3qTHfLII9DYCM88k3UkZtZZNDY20tjY2OHvadfVYpLOjohrOrz15Ht3B2ZFxP5poiEiLk/XzSY5\n5fUSMDci+qXtxwNDImJM2qchIuZJ6gKsiIhdWthOVV0ttn49DBoE558PI0dmHY2ZdVbleJ7L6pzL\nkMe3dUO5JNXlvB0OLEyXZwKjJG0jaU+gDpgfEa8CayQNTgv8JwL35HymaSR1HDCnI7FViptugo98\nBL7ylawjMTNru7ZOuV8sl0naB9gA/AUYAxARSyTdDiwhuft/bM5wYyzJXGbdSEY6s9P2m4FpkpYC\nK0nmPqtq//M/8P3vJ6fEXMQ3s2rUrin3Jb0cEX3z96wc1XRa7NRTYccd4aqrso7EzDq7ok+5L+kf\nJFPtt8S38pXIvHlw//0u4ptZddticomI7be0zkpjwwY4/XSYNAk+/OGsozEza7+CC/pWev/+78n0\nLieckHUkZmYd48ccV4iVK6F//+QpkwMGZB2NmVmivTUXJ5cK8e1vwzbbwHXXZR2JmdlmRS/oW/n8\n8Y9wzz0u4ptZ7XDNJWMbNybT6V96KeywQ9bRmJkVh5NLxm69NblR8qS8s7WZmVUP11wy9MYb0K8f\n3HdfMo+YmVmlcUE/j0pMLmeemUxQeeONWUdiZtYyF/SrzKJFcPvtsGRJ1pGYmRWfay4ZiEiK+Bdf\nnMx8bGZWa5xcMvCLX8C77yZPmDQzq0WuuZTZ6tVJEf+uu2Dw4KyjMTNrnQv6eVRKcvmXf4E1a5J5\nxMzMKp0L+lXgT39KTok9/XTWkZiZlZZrLmUSkVx6fNFFsMsuWUdjZlZaTi5lMmNGctPkt76VdSRm\nZqXnmksZvPlmUsSfMQO+8IVMQjAza5f21lwyGblIuljSk5IWSZojqW/avoekdyQtTF835HxmkKTF\nkpZKujanvaukGWn7PEkfz2KfWvPDH8LQoU4sZtZ5ZDJykdQjIt5Ml88EBkTENyTtAdwbEQe08Jn5\nwBkRMV/SLOC6iJgtaSywf0SMlTQSOCYiRrXw+UxGLs8+C4ceCosXQ69eZd+8mVmHVNXIpSmxpLYH\nXmutv6TeQI+ImJ82TQVGpMvDgCnp8h3A0CKG2iFNRfyJE51YzKxzyaygL+kSSX8FTgIuz1m1Z3pK\nrFHSIWlbH2BZTp/laVvTupcBImI9sFrSTqWNvjB33gkrVsDpp2cdiZlZeZXsPhdJDwIt/b1+fkTc\nGxETgYmSzgMmA6cArwB9I+INSQcBd0var1gxNTQ0bFqur6+nvr6+WF/9AW+/ndwwOWUKbL11yTZj\nZlZUjY2NNDY2dvh7Mr9aTNLuwKyI2L+FdXOB8cAK4OGI6Je2Hw8MiYgxkmYDDRExT1IXYEVEfOBO\nknLXXC64AP7yF5g+vWybNDMruqqquUiqy3k7HFiYtu8saat0eS+gDng+IlYAayQNliTgROCe9PMz\nSU6tARwHzCnDLrRq6VK46Sa48sqsIzEzy0ZW079cJmkfYAPwF2BM2j4E+IGkdcBG4FsRsSpdNxa4\nFehGMtKZnbbfDEyTtBRYCXzgSrFyioCzzoJzz4U+ffL3NzOrRZmfFiuXcp0WmzkzSSxPPgnbbFPy\nzZmZlZQnrqwA77wDZ58NP/mJE4uZdW6eW6yIfvQjOOggOPzwrCMxM8uWT4sVyQsvwMEHw4IFsPvu\nJduMmVlZVdXVYrXoO99J7mtxYjEzc82lKH7zm+QBYDNmZB2JmVll8Milg957D8aNg+uug65ds47G\nzKwyOLl00FVXQf/+cNRRWUdiZlY5XNDvgL/+FQYOhD/+Efbcs6hfbWZWEVzQz8D48cmU+k4sZmbv\n54J+Oz30EDzxBEydmnUkZmaVxyOXdli7NhmxXHMNdOuWdTRmZpXHyaUdrr02ORV29NFZR2JmVplc\n0G+j5cthwAB47DGoq8vf38ysmrW3oO/k0kZf/WoyarnkkiIEZWZW4Zxc8ihGcvntb+HrX4clS2C7\n7YoUmJlZBfOlyCW2bh2ccUZy06QTi5lZ65xcCvTjH0OvXnDssVlHYmZW+XxarACvvgr77w+PPgr7\n7lvkwMzMKlhVnhaTNF7SRkk75bRNkLRU0rOSjshpHyRpcbru2pz2rpJmpO3zJH282HGeey6ceqoT\ni5lZoTK7Q19SX+Bw4KWctv7ASKA/0Ad4SFJdOuS4ERgdEfMlzZJ0ZETMBkYDKyOiTtJIYBIwqlhx\n/v73MGcOPPNMsb7RzKz2ZTlyuRr4XrO24cD0iFgXES8CzwGDJfUGekTE/LTfVGBEujwMmJIu3wEM\nLVaAGzYkRfwrroAePYr1rWZmtS+T5CJpOLAsIp5qtmpXYFnO+2UkI5jm7cvTdtKfLwNExHpgde5p\nto74yU+gZ08YVbRxkJlZ51Cy02KSHgR6tbBqIjABOCK3e6niaK+//x0aGuDhh0EVF52ZWWUrWXKJ\niMNbape0P7An8KSS39q7AU9IGkwyIumb0303khHL8nS5eTvput2BVyR1AXpGxOstbbuhoWHTcn19\nPfX19VuMf8IEOOGE5CoxM7POorGxkcbGxg5/T+aXIkt6ARgUEa+nBf3bgM+QFvSBvSMiJD0OjAPm\nA/cB10XEbEljgQMiYoykUcCIiPjAiay2XIo8fz6MGJEU8Xv2LMpumplVpfZeilwJz3PZ9Bs/IpZI\nuh1YAqwHxuZkhLHArUA3YFZ6pRjAzcA0SUuBlXTwSrENG+D00+Hyy51YzMzaK/ORS7kUOnL52c9g\nyhR45BHXWszMPHFlHoUkl9dfh3794P774cADyxSYmVkFc3LJo5DkMmYMbLUVXH99mYIyM6tw1Vxz\nqQgLFsBdd/lOfDOzYvCsyMDGjUkR/5JLYMcds47GzKz6ObkAU6cmCeaUU7KOxMysNnT6msuqVUkR\n/9574dOfziAwM7MK5oJ+HltKLmedBe++m8wjZmZm7+eCfjs89RRMn+4ivplZsXXamktEUsT/wQ/g\nIx/JOhozs9rSaZPLbbfB22/DaadlHYmZWe3plDWXNWuSIv5//id87nMZB2ZmVsFc0M8jN7mcc04y\n1cstt2QclJlZhXNyyaMpuSxZAocdBk8/DR/9aNZRmZlVtvYml05Vc4mAM8+ECy90YjEzK6VOlVx+\n9St47bVkgkozMyudTnVabLfdgttug0MPzToaM7Pq4NNiBaivd2IxMyuHTjVyeeWVoHfvrCMxM6se\nvlosj0Ifc2xmZptV5WkxSeMlbZS0U/p+D0nvSFqYvm7I6TtI0mJJSyVdm9PeVdKMtH2epI9nsS9m\nZrZZZslFUl/gcOClZquei4iB6WtsTvuNwOiIqAPqJB2Zto8GVqbtk4FJpY69EjU2NmYdQkl5/6pX\nLe8b1P7+tVeWI5erge8V0lFSb6BHRMxPm6YCI9LlYcCUdPkOYGgxg6wWtf4/uPevetXyvkHt7197\nZZJcJA0HlkXEUy2s3jM9JdYo6ZC0rQ+wLKfP8rStad3LABGxHljddJrNzMyyUbLnuUh6EOjVwqqJ\nwATgiNzu6c9XgL4R8Yakg4C7Je1XqhjNzKw0yn61mKT9gTnA22nTbiQjkc9ExP806zsXGA+sAB6O\niH5p+/HAkIgYI2k20BAR8yR1AVZExC4tbNeXipmZtUNVPIkyIv4EfKzpvaQXgEER8bqknYE3ImKD\npL2AOuCcpacqAAAHCUlEQVT5iFglaY2kwcB84ETguvQrZgInAfOA40gSV0vbbfN/HDMza59KeMxx\n7ohiCPADSeuAjcC3ImJVum4scCvQDZgVEbPT9puBaZKWAiuBUWWJ2szMtqjT3ERpZmblU3Nzi0k6\nUtKz6U2V526hz3Xp+iclDSx3jB2Rb/8k1UtanXMj6gVZxNkekm6R9DdJi1vpU83HrtX9q/Jj11fS\nXElPS/qTpHFb6FeVx6+Q/avy47etpMclLZK0RNJlW+hX+PGLiJp5AVsBzwF7AFsDi4B+zfp8keS0\nGsBgYF7WcRd5/+qBmVnH2s79OxQYCCzewvqqPXYF7l81H7tewIHp8vbAf9fYv71C9q9qj18af/f0\nZxeSGvYhHTl+tTZy+QzJHf4vRsQ64D+A4c36bLrpMiIeB3aQ9DGqQyH7B5sv7a4qEfEI8EYrXar5\n2BWyf1C9x+7ViFiULv8DeAbYtVm3qj1+Be4fVOnxA4iIpit4tyH5Q/b1Zl3adPxqLblsuqEytYzN\nN1u21me3EsdVLIXsXwCfT4etsyT1L1t0pVfNx64QNXHsJO1BMkJ7vNmqmjh+rexfVR8/SR+StAj4\nGzA3IpY069Km41cJV4sVU6FXJzT/66JarmooJM4FJDeivi3pKOBu4JOlDausqvXYFaLqj52k7YH/\nBM5K/8L/QJdm76vq+OXZv6o+fhGxEThQUk/gfkn1EdHYrFvBx6/WRi7Lgb457/vy/mljWurTdBNn\nNci7fxHxZtPwNiJ+A2xdQ9PhVPOxy6vaj52krUnm9/tFRNzdQpeqPn759q/aj1+TiFgN3Ad8utmq\nNh2/WksufySZMXkPSdsAI0lussw1E/g6gKTPAqsi4m/lDbPd8u6fpI9JUrr8GZLLzZufO61W1Xzs\n8qrmY5fGfTOwJCKu2UK3qj1+hexflR+/nSXtkC53I5mxfmGzbm06fjV1Wiwi1ks6A7ifpCB1c0Q8\nI+lb6fqfRMQsSV+U9BzwFnBKhiG3SSH7RzJLwRhJ60mm2Kmam0olTQcOA3aW9DJwEclVcVV/7CD/\n/lHFxw74AvA14ClJTb+Uzgd2h5o4fnn3j+o+fr2BKZI+RDLomBYRczryu9M3UZqZWdHV2mkxMzOr\nAE4uZmZWdE4uZmZWdE4uZmZWdE4uZmZWdE4uZmZWdE4uZu0kqVHSoDJsZ1w6Dfq0Em7jxWq8m9wq\nV03dRGlWZu2+SUxSl4hYX2D3McDQiHilvdsrgG94s6LyyMVqWjpVzjOSfpo+5Ol+Sdum6zaNPNLp\nL15Il0+WdLekByS9IOkMSedIWiDpMUk75mzixPTBUIslHZx+fjslDwZ7PP3MsJzvnSlpDvBgC7H+\nS/o9iyWdlbbdBOwFzJZ0drP+J0u6R8lDrP4s6cKcdV9Lt79Q0k3pnddIukHSH9L/Fg0txNBN0m8k\njZbUXdJ9Sh4gtVjSVzpwKKyTcXKxzmBv4PqI2B9YBRybtgdb/ot9P+AY4GDgEmBNRBwEPEY6vxLJ\nDLHdImIgMBa4JW2fCMyJiMHA/wGukNQ9XTcQODYi/il3Y2mSO5nkmT2fBU6TNCAivg28AtRvYU6r\ng4EvAZ8CvixpkKR+wFeAz6exbQROaIotIg4GBgCHSdo/57t6kMwf9cuIuBk4ClgeEQdGxAHA7C38\ntzL7AJ8Ws87ghYh4Kl1+guRJnvnMjYi3gLckrQLuTdsXk/wihyQxTYfkQWCSPqxkuvIjgKMlnZP2\n60oyB1UAD0bEqha2dwhwZ0S8AyDpTmAI8GSeOB+IiDdyPnMIsAEYBPwxnUexG/Bq2n+kpNNI/u33\nBvoDfyJJlPcAkyJietr3KeBKSZcDv46IR/PEYraJk4t1Bu/lLG8Atk2X17N59L4t75f7mY057zfS\n+r+bppHQlyJiae4KSYNJJvzb0udyn5Uh8tdBmq/P/cyUiDi/2fb3BMYDn46I1ZJ+zub9DuBRktFK\nU8JcquQ56f8M/FDSnIi4OE9MZoBPi1nn1PRL/EU2P7PiuDZ+tml5JICkQ0imIF9DMmv1uE2dkl/Q\nzT/b3CPAiLTmsR0wIm3LF8vhknZMp0kfTpIg5gDHSdol3f5OknYnOe31FrBGyeNpj2r2fRcCb0j6\ncfq53sC7EfFL4ErgoDzxmG3ikYt1Bs3/wm96fyVwu6RvkjwcKXLWRwv9m68L4F1JC0j+LZ2atl8M\nXCPpKZI/4J4nef74Fms8EbFQ0q3A/LTpZxHRdEpsSyOYSPvfQfLgpmkRsQBA0gXAA2khfx0wNiLm\np9PFP0vyuNoPnOaKiLPSixEmkSSpKyRtBNaSXLVmVhBPuW9WpSSdDAyKiDOzjsWsOZ8WM6terV3t\nZpYpj1zMzKzoPHIxM7Oic3IxM7Oic3IxM7Oic3IxM7Oic3IxM7Oic3IxM7Oi+18U50dFyFSy7AAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1095c0fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(log_evidence[-4:])\n",
    "plt.ylabel('Log-evidence')\n",
    "plt.xlabel('number of peaks')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Bayes factor is related to the exponential of the difference between the log-evidence values.  Thus, 0 peaks is not very likely compared to 1 peak. But 1 peak is not as good as 2 peaks. 3 peaks is not that much better than 2 peaks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r01 = np.exp(log_evidence[-4] - log_evidence[-3])\n",
    "r12 = np.exp(log_evidence[-3] - log_evidence[-2])\n",
    "r23 = np.exp(log_evidence[-2] - log_evidence[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 1.43158716924e-50 0.909965654827\n"
     ]
    }
   ],
   "source": [
    "print(r01, r12, r23)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These numbers tell us that zero peaks is 0 times as likely as one peak. Two peaks is 7e49 times more likely than one peak. Three peaks is 1.1 times more likely than two peaks. With this data one would say that two peaks is sufficient. Caution has to be taken with these values. The log-priors for this sampling are uniform but improper, i.e. they are not normalised properly.  Internally the lnprior probability is calculated as 0 if all parameters are within their bounds and `-np.inf` if any parameter is outside the bounds. The `lnprob` function defined above is the log-likelihood alone. Remember, that the log-posterior probability is equal to the sum of the log-prior and log-likelihood probabilities. Extra terms can be added to the lnprob function to calculate the normalised log-probability.  These terms would look something like:\n",
    "\n",
    "$\\log \\left(\\prod_i \\frac{1}{\\max_i - \\min_i}\\right)$\n",
    "\n",
    "where $\\max_i$ and $\\min_i$ are the upper and lower bounds for the parameter, and the prior is a uniform distribution. Other types of prior are possible. For example, you might expect the prior to be Gaussian."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
